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I.  Introduction 


This  is  a  report  on  the  feasibility  of  determining  atmospheric  conditions  from  a  pro¬ 
jectile’s  observed  trajectory.  In  this  preliminary  study,  ‘measurements’  are  generated  as 
solutions  of  the  modified  point-mass  (MPM)  equations  of  motion,  using  the  1962  U.S.  Stan¬ 
dard  Atmosphere  and  specified  wind  profiles.  The  hope  is  that  eventually  the  generated 
data  will  be  replaced  in  the  program  by  instrument  readings  of  trajectory  variables. 

The  over-all  scheme  is  as  follows.  The  MPM  equations  are  altered  by  re-defining  the 
air  density,  sound  speed  and  wind  in  terms  of  a  small  number  of  parameters.  A  nonlinear 
curve-fitting  technique  known  as  FINLIE  (Ref.  1]  then  fits  the  revised  differential  equations 
to  the  generated  data  by  adjusting  the  parameter  values.  The  original  and  fitted  air  density, 
sound  speed  and  wind  histories  can  then  be  compared. 

The  FINLIE  user  must  construct  a  FORTRAN  subroutine  defining  the  original  set 
of  equations  plus  an  auxiliary  set.  The  auxiliary  equations  involve  partial  derivatives  of 
the  dependent  variables  of  the  original  set.  Expressions  for  these  partial  derivatives  can 
be  extremely  complicated  and  the  opportunity  for  error  -  both  in  deriving  the  expressions 
and  in  coding  them  -  is  great.  Hence  a  software  package  called  MACSYMA'  [Refs.  2  and 
3]  was  used.  From  the  original  equations,  this  versatile  program  determined  expressions 
for  each  of  the  partial  derivatives  and  presented  these  expression  as  FORTRAN  code. 

In  the  next  few  sections,  we  will  discuss  the  modified  point-mass  equations,  the  stan¬ 
dard  atmosphere  used,  the  revised  equations  used  to  fit  the  data,  and  the  application  of 
FINLIE/MACSYMA  to  the  fitting  problem.  Succeeding  sections  will  discuss  results  and 
conclusions. 


II.  The  Modified  Point-Mass  Equations 


The  modified  point-mass  trajectory  model  is  the  primary  method  of  simulating  tra¬ 
jectories  in  the  preparation  of  Firing  Tables.  The  basic  equation  of  motion,  derived  in  Ref. 
4.  can  be  re-written  [Ref.  5]  in  the  form 

U  ^  (D  -  A)V  +  Ga\  (1) 

where 


D 

A 

G.x 


ha(G*V) 

(1  +/la)F2 


*  MACSYMA  is  a  trademark  of  Synibolie.s,  Iin  ..  '257  V'assar  St.,  Cambridge.  MA  ()'213y 
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and  where 

U  =  projectile  velocity  with  respect  to  the  earth 
V  =  U  —  W  =  projectile  velocity  with  respect  to  the  air 
\y  =  wind  velocity  with  respect  to  the  earth 
G  =  the  sum  of  the  gravity  and  Coriolis  accelerations 
(p  =  axial  spin,  rad/s. 

(All  symbols  are  defined  in  the  List  of  Symbols.) 

The  ajcial  spin  (t>  is  obtained  as  the  solution  of  a  simplified  roll  equation: 


where 


(2) 


The  aerodynamic  coefficients  and  Cg^,  in  Eqs.(  1)  and  (2)  are  assumed 

to  be  functions  of  Mach  number.  The  drag  coefficient  Co  is  assumed  to  depend  on  both 
Mach  number  and  the  yaw  of  repose,  cTe.  In  particular,  Cp  is  assumed  to  have  the  form 


Cp  —  Cpo  +  Cp2\Qe 


where  Cpo  is  a  function  of  Mach  number  and  Cp2  is  a  constant. 

The  yaw  of  repose  can  be  computed  from  the  relation  [Ref.  5]: 


cTe 


^GaxV 
BV^Cm^  • 


(3) 


(4) 


III.  Component  Form 


The  vector  equation  (1)  is  not  in  a  suitable  form  for  programming;  we  need  component 
expressions.  Hence  our  first  task  is  to  select  a  convenient  coordinate  system  for  describing 
the  motion  of  a  projectile  along  its  trajectory. 
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For  the  moment,  assume  that  the  launch  point  is  at  some  sea-level  spot  on  the  earth’s 
surface.  Then  we  set  our  origin  at  the  launch  point  and  define  a  right-handed  Cartesian 
system  as  follows:  the  1-  and  3-axes  form  a  plane  tangent  to  the  earth  at  the  origin;  the 
2-axis  is  perpendicular  to  this  plane,  positive  upward,  and  the  1-axis  is  chosen  so  that 
the  velocity  U  at  time  zero  is  in  the  1-2  plane.  We  will  call  this  coordinate  system  the 
‘flat-earth’  system.  Then  the  projectile’s  position  vector  X  with  respect  to  the  earth  can 
be  written  in  component  form  as 

X  =  (X,,  .Y2,  A'3) 

where 


A'l  is  the  down-range  distance, 

A'2  is  the  height  above  the  1-3  plane, 

A’3  is  the  lateral  distance,  positive  to  the  right  when  looking  down-range. 

and  where  A"  =  U .  Since  the  trajectory  -  in  most  cases  -  lies  nearly  in  the  T2  plane.  A*3 
is  usually  much  smaller  in  magnitude  than  A'l  and  A'2. 

Similarly,  we  can  write 

U  =  iUu  u,,  ih) 

V  =  (V'l,  V2,  F3) 

w  =  (W,,  W2,  W3) 

G  =  (Gi,  G2,  G3). 

where  ['3  is  usually  much  smaller  in  magnitude  than  U\  and  ^2-  The  initial  velocity  is 
given  by 

Uo  =  |L',5|(cos £■,  sinE,  0) 

where  E  is  the  gun  elevation. 

In  the  system  we  have  just  described,  the  launch  point  is  at  (0.0.0).  However,  if 
the  launch  point  is  not  at  sea  level  (for  example,  if  it  is  half-way  up  a  mountain),  then 
we  choose  our  origin  at  the  sea-level  point  directly  beneath  the  launch  point.  Thus  the 
origin  of  our  system  may  be  under  the  earth’s  surface.  In  general,  the  launch  point  is  at 
(0.  .Y90,  0),  where  A'20  is  the  height  of  the  launch  point  above  sea  level.  The  motivation 
for  this  rather  cumbersome  shifting  of  the  origin  is  to  retain  A'2  as  a  measure  of  the  height 
above  sea  level  for  a  flat-earth  model. 


Eq.(l)  can  be  written  in  component  form  as 


(-2 

f'3 


{D-A)\\  -k 


(D-A)\2  + 


1 

\  +  ha 
1 

1  ha 


{ D  —  .4 )  V  3  -l- 


1 

1  -f-  ha 


G. -f 
G2  -f 
G  3  4- 


hL{  G2\'j  —  G3V2  ) 
il-hM)V 
hL(G,l\-G,V3) 

( 1  -  h,\i  )\ 
ftLiGilj  —  G2I  1  ) 
{l-hM)V 


(0) 

(6) 


.1 


Before  we  can  program  these  equations,  we  need  expressions  for  the  components  of 
G.  Recall  that  G  is  the  sum  of  the  gravity  and  Coriolis  accelerations: 

G  =  g  +  C  (8) 

The  gravity  vector  can  be  written  as 

go{Xr/R,  l  +  X^/R,  X3/R) 

^  [(A'i/R)2  +  (I  +  .Y2W  +  (.Y3/i?)2p/2  ^  ^ 

where 

go  =  magnitude  of  g  at  the  origin  (the  standard  value  at  45  deg.  N.  latitude  being 
9.S0665  m/s^) 

R  =  'effective'  radius  of  the  earth  (6  356  766  m). 

This  expression  for  g  is  usually  .simplified  by  ignoring  X3/R  and  assuming  that  Xi/R 
and  X2/R  -  while  not  ignorable  -  are  much  less  than  unity.  Then  we  have 

^  go{X\/R,  l+Xj/R,  0) 

^  ~  ~  (l+A2/i?P 

«  -go{Xt/R,  {1  +  X2/R)~\  0) 

^  -go(Ai/R,  I-2.Y2/R,  0)  (10) 

The  Coriolis  acceleration  is  given  by 

C  =  ~2wxU  (11) 

where  w  is  the  earth's  angular  velocity  vector: 

u;  =  (wi,  u;2,  0:3) 

=  cj£;(cosT  COS-4Z,  sin£,  —cosL  sinTZ)  (1*2) 


where 


=  |w| 


27r  rad/ sidereal  day 
86164.09  s/sidereal  day 


ss  7.291  X  10  ^  rad/s 


and  where 

L  =  latitude  at  launch  point  (for  the  Southern  Hemisphere,  replace  L  by  -L). 
AZ  =  azimuth  of  the  l-axis,  measured  clockwise  from  North. 

Hence  we  have 


Cl  —  —goXi/R  —  2(u;2£3  —  u-’it  2) 

C2  =  -^7„(l  -  2X2/ R)  -  2(^’3C,  -  ^•1^3) 

C3  =  — 2(j;i £2  ~  ^'2^  1 ) 


•t 


(13) 

(14) 

(15) 


IV.  Atmospheric  Conditions 


The  equations  of  motion  used  in  this  investigation  are  essentially  Eqs.(5-7)  and  (2). 
However,  two  versions  of  this  set  of  equations  had  to  be  programmed:  one  to  generate  data 
(A'l.  .\'2,  and  .Y3  vs.  t)  and  one  to  fit  the  generated  data. 

In  the  generating  version,  the  wind,  air  density,  and  speed  of  sound  are  somewhat 
complicated  functions  of  altitude.  We  will  discuss  these  functions  shortly. 

In  the  fitting  version,  simpler  expressions  are  assumed  for  the  wind,  air  density,  and 
sound  speed.  We  will  see  later  that  these  expressions  involve  just  five  fitting  parameters. 
The  whole  point  of  this  exercise  is  to  fit  the  second  version  of  the  equations  to  the  given 
trajectory  (probably  piecewise)  to  obtain  values  for  those  five  parameters. 


1.  Atmosphere  Used  in  the  Generating  Equations 

For  purposes  of  generating  test  trajectories,  we  used  the  air  density  and  speed  of  sound 
from  the  U.S.  Standard  Atmosphere,  1962  [Ref.  6],  hereafter  called  STAT.  The  ‘geometric’ 
altitude  Z  discussed  in  STAT  is  not  quite  the  same  as  .Y2,  the  height  above  the  1-3  plane 
in  our  flat-earth  system.  A  better  approximation  to  the  altitude  allows  for  the  curvature 
of  the  earth.  We  have^ 

(R  +  X2f  +  X^,={R  +  Zf 

so  that 


Z  =  ,JiR  +  X2)^  +  X^-R 
^  X2  +  X^J2R 


:i6) 


Note  that  in  approximation  (16),  the  difference  between  Z  and  A'’2  is  solely  a  function  of 
the  down-range  distance  A'p  At  A'l  =  25  kilometers,  for  example.  Z  —  A'l  is  about  50 
metres  for  any  value  of  A’2.  If  a  trajectory  is  to  be  run  to  impact  at  ground  level,  we 
should  set  A"2  =  —Xi/2R,  not  A’’2  =  0,  as  the  stopping  point. 


The  STAT  value  for  g.  the  magnitude  ot  g,  is  given  by 

9o 


9i^)  = 


{1+Z/Rr 


(IT) 


In  the  STAT  system,  the  basic  altitude  (up  to  Z  =  90  km)  is  the  'geopotential'  altitude 
H,  defined  as 

g(s)  ds  Z 


H 


-i: 


IS) 


go  I  +  Z/R 

ST.4T  divides  altitude  H  into  eight  zones,  in  each  of  which  the  temperature  is  assumed  to 


be  li 


near 


in  H: 


T  =  TB  +  r{H  -  Hh) 


{deg  K 


19) 


^  .‘"'trictly,  th#*  geometric  altitutle  Z  i.s  ineasureci  /iloiig  a  line  of  gravitational  force:  hf»wever,  the  distinction  hetwef'n  this  line 
)f  foreo  and  a  straight  line  Is  negligible  for  our  pur|>oses. 


Table  1.  1962  US  Standard  Atmosphere  Values 


Zone 

H  (km) 

Hb  (km) 

Tb  (deg  K) 

T'  (deg  K/km) 

Pb  (kg/m-’) 

1 

-5,  11 

0 

288.15 

-6.5 

1.225000+0 

2 

11,  20 

11 

216.65 

0 

3.6391803-1 

3 

20,  32 

20 

216.65 

1 

8.8034864-2 

4 

32,  47 

32 

228.65 

2.8 

1.3225009-2 

5 

47,  52 

47 

270.65 

0 

1.4275335-3 

6 

52,  61 

52 

270.65 

-2 

7.5943220-4 

7 

61,  79 

61 

252.65 

-4 

2.5109063-4 

8 

79,  89 

79 

180.65 

0 

2.0011372-5 

where  Tb  is  the  value  of  T  at  Hb,  the  bottom  altitude  for  the  zone  (except  in  the  lowest 
zone’,  where  Hb  =  0  but  the  zone  extends  down  to  -5  kilometres).  Values  of  Hb-  Tb-  and 
the  slope  T'  (  =  dT/dH)  are  given  in  Table  1  for  each  of  the  eight  zones. 

The  speed  of  sound  at  a  given  altitude  can  be  computed  from  the  corresponding 
temperature  (deg  K); 

V;  =  20.04680276  \/T  {m/s)  (20) 


For  those  zones  in  which  T'  is  not  zero,  the  STAT  air  density  is  given  by 


£_ 

PB 


j,  -i-Q/r 

\f^) 


(21) 


where 


Q- 


go  *  molecular  weight 


34.16319474  (deg  K/km) 


universal  gas  constant 
and  where,  as  before,  subscript  B  denotes  the  value  at  H  =  Hb- 

For  those  zones  in  which  T'  =  0  (that  is.  in  which  the  temperature  is  assumed  con¬ 
stant),  the  STAT  air  density  is  given  by 


p!  Pi 


-Q{H~HB)/Tg 


{22\ 


We  obtained  the  approximate  values  of  ps  listed  in  Table  1  by  setting  ps  =  1.225 
kg/m^  for  zone  1  and  then  calculating  pb  for  each  succeeding  zone  as  the  density  at  the 
top  of  the  previous  zone. 

There  is  no  'standard'  wind  profile;  hence  the  winds  used  in  generating  trajectories 
were  arbitrarily  chosen  as  follows:  the  vertical  wind  component.  VV’2'  was  zero  while  TUi 
and  IV’'5  were  specified  quadratics  in  height  A*2. 

*  “Weather”,  as  wp  usually  think  of  it  (clouds,  rain,  snow  etc.)  is  almost  always  confined  -  like  most  of  us  -  to  zone  1.  the 
troposphere.  Although  the  nominal  rate  of  temperature  flerrease  in  zone  1  is  H..5  deg  K/krn,  loiral  temperature  inversions  are 
f'om  rn(')n. 
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2.  Atmosphere  Used  in  the  Fitting  Equations 


Eventually  (but  not  in  this  preliminary  study)  we  intend  to  divide  the  traif  ctory  into 
overlapping  time  inter\als,  some  arbitrary  time  within  one  interval  (say,  th  .  mid-time) 
serving  as  the  initial  time  for  the  next  segment.  (These  time  intervals  are  not  to  be 
confused  with  STAT  altitude  zones;  a  time  interval  could  lie  entirely  in  one  altitud  )ne 
or  involve  two  or  more  zones.) 

For  the  fitting  equations,  we  decided  that  over  any  segment  of  the  trajectory  the  wind 
could  be  approximated  by  the  relations 

\V\  =  Cx+C2X2 
W2  =  0 
IVa  =  C3  +  C,X2 


(23) 

(24) 


This  gave  us  four  parameters  (so  far)  to  be  determined  by  the  fit:  C]  through  C.i.  To  keep 
the  total  number  of  parameters  to  a  minimum  (for  reasons  that  will  become  clearer  when 
we  discuss  FINLIE),  we  decided  to  express  the  speed  of  sound  and  the  air  density  in  terms 
of  a  single  parameter: 

C5=T'  (26) 

The  speed  of  sound  is  assumed  to  have  the  form  of  Eq.(20): 


U3  =  20.04680276  ^Ta  +  C^iH  -  Ha) 


(m/s) 


(27) 


where  Ta  and  Ha  -  the  temperature  and  altitude  at  the  starting  time  of  the  interval  under 
consideration  -  are  assumed  to  be  known  quantities. 

The  density  relation  (21)  is  unsuitable  for  fitting  purposes  because  it  tends  to  become 
indeterminate  as  T'  (=  C5)  goes  to  zero.  Actually,  in  an  earlier  phase  of  this  investigation, 
we  did  use  Eq.(21)  to  fit  a  generated  trajectory  that  remained  in  zone  1;  that  is,  for  which 
the  altitude  was  less  than  11  km.  Starting  from  a  non-zero  estimate  for  C5,  FINLIE  zeroed 
in  on  the  correct  value  (-6.5  deg  K/km)  in  a  few  iterations.  However,  we  can't  in  general 
be  sure  that  C.5  will  be  nonzero  in  the  fitting  interval  and  we  certainly  don't  want  to 
preclude  an  initial  estimate  of  zero.  Hence  we  decided  to  re-write  Eq.(21)  so  as  to  remove 
the  indeterminacy.  To  do  so.  we  set 

^  TB-kTjH-HB)  ^  ,  V 

Tb  Tb 

where 

^  _  T{H-Hb) 

Tb 

We  s('e  from  Table  1  that  in  any  zone.  |A|  is  less  than  1;  hence  we  can  express  ln(  1  -t-  A) 
as  a  convergent  power  series: 

A'^  A"’  A' 

\n{T/TB)  =  A~~  +  ~-  —  +  ---. 


Eq.(21)  can  then  be  written  as 
P  , 


=  exp[-{l+QIT)\n{T/TB)\ 

PB 

f  A  A^ 

=  exp  -(l  +  g/r)A(l--  +  — -  — +  ■ 


exp  — 


{r  +  Q){H-HB) 


A  A^  A^ 


This  form  of  the  STAT  density  expression  (21)  has  the  computational  advantage  that  it 
offers  no  difficulties  as  T'  goes  to  zero;  in  fact,  at  T'  =  0,  the  new  form  reduces  to  Eq.(22). 
In  our  fitting  equations,  we  will  eventually  use  a  truncated  form  of  Eq.(2S): 

pIpa  =  exp[-{C5  +  Q)iH  -  Ha)/Ta\  (29) 

where  pa  -  the  value  of  p  at  the  start  of  the  fitting  interval  -  is  assumed  to  be  known.  For 
this  preliminary  study,  we  used  the  expanded  form  (28)  in  order  to  recover  -  if  possible  - 
the  generated  value  of  T' . 

Only  the  first  time-interval  values  of  Ta,  Ha,  and  pa  (presumably  obtained  in  ‘real  life' 
by  measurements  at  the  launch  site)  ar<“  required  inputs  to  the  fitting  process:  thereafter, 
the  closing  values  for  the  k-th  interval  become  the  starting  values  for  the  (k+l)-th  interval. 


V.  FINLIE 

FINLIE  [Ref.  1]  is  a  FORTRAN  program  for  fitting  a  system  of  first-order  ordinary 
differential  equations  to  measurements  taken  on  one  or  mon  of  the  dependent  variables. 
No  knowledge  of  the  system’s  solution  is  required;  indeed,  in  most  cases,  no  closed-form 
solution  exists. 

Our  first  step  in  using  FINLIE  was  to  re-write  the  equations  in  the  required  form.  We 


y\  =  -W 

;/4 

=  -V,  =  Fi 

y2  =  -^"^2 

yh 

II 

II 

1/3  =  -Y3 

ye 

—  -Y3  —  ^'3 

yr 

=  0 

7)  and  (2)  could  be 

e.xpressed  as  a 

system  of  seven  first-order  equations; 

y\  =  E,  [=  y4] 

m  = 

F,  [=  RHS  of  Eq.(o)\ 

ih  =  Fo  [=  ys] 

y's  = 

F5  [^RHSof  Eq.[^)] 

y-i  =  Fi  [=  yr,] 

ye  = 

Fe  [=RHS  of  Eq.(-)\ 

(30) 

yr  = 

F,  [=  -BC,^[y:CIV)] 

Tlie  FINLIE  user  must  construct  a  FORTRAN  subroutine  defining  the  system  of 
equations  -  in  our  ca,se.  system  (30)  -  and  an  auxiliary  set  of  so-called  influence  equations. 
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To  describe  the  latter  set,  we  introduce  qk  where 


<i\,  ■■■■,  Q7  —  yio?  •••1  yjo 
98)  •■•)  9i2  =  Cl, Cs 

Influence  coefficients,  the  dependent  variables  of  the  influence  equations,  are  defined  as 

(31) 

dqk 

where  index  j  runs  from  1  to  n  and  index  k  runs  from  1  to  (n  +  N),  n  being  the  order  of 
the  original  system  (7  in  our  case)  and  N  the  number  of  parameters  (here  5).  Then  the 
influence  equations  can  be  written  in  the  form  [Ref.  1]: 

=  E  (^1  D.,  +  C,,  (32) 


where 


Cjk  =  0 

dF, 

9Ck-n 


if  k  <  n 
if  k  >  n 


FINLIE  automatically  assigns  the  proper  initial  conditions  to  the  influence  equations 
[(Djk)o  is  1  if  j  =  k  and  0  otherwise]  and  integrates  the  influence  equations  simultane¬ 
ously  with  the  original  system. 

Note  that  there  are  n(n  +  N)  influence  equations.  This  can  be  an  uncomfortably  large 
number.  We  can’t  usually  do  anything  about  n  but  it  helps  if  we  can  keep  N.  the  number 
of  parameters,  as  small  as  possible.  This  is  why  we  restricted  ourselves  to  5  parameters 
and  'only’  7(7  +  5)  =  84  influence  equations. 

Because  of  the  simplicity  of  the  first  three  of  the  seven  equations  in  set  (30),  the 
corresponding  influence  equations  are  equally  simple: 

Dxk  =  D,k  (33) 

D2k  =  D^k  (34) 

Dsk  =  Ddk  (35) 

The  remaining  four/seventh  of  set  (32)  is  not  quite  so  easily  expressed.  Consider,  for 
example,  just  one  of  the  partial  derivatives  needed  to  evaluate  D^k,  say 


where  is  the  right-hand  side  of  Eq.(5).  Note  that  A’^  enters  F.i  in  three  ways: 

•  through  the  gravity  component  G2,  Eq.(14); 

•  through  the  geometric  altitude  Z,  Eq.(16).  and  from  there  into  H  and  p  and 
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•  through  the  winds.  ITi  and  IV3  [Eqs.(23)  and  (25)],  and  from  there  into  velocity  V. 


And,  of  course,  since  the  aerodynamic  coefficients  are  specified  functions  of  Mach  number 
(=  V/V,),  A'2  is  involved  wherever  there  is  an  aerodynamic  coefficient.  Thus  determining 
an  expression  for  dFildX2  -  or  indeed  for  many  of  the  other  partials  -  is  not  a  simple 
task. 


Two  factors  came  to  our  aid  at  this  point.  First,  we  were  able  to  make  some  simplifying 
assumptions.  As  Ref.  1  points  out, 

...  certain  liberties  can  be  taken  with  the  influence  equations:  expressions  can  be 
approximated  by  simpler  ones,  the  effect  of  certain  [parameters]  on  certain  terms 
in  the  original  equations  can  be  ignored,  etc.  If  done  with  care  and  judgment, 
such  simplifications  will  have  no  effect  on  the  final  answer:  the  same  [fit]  will  be 
reached  with  or  without  the  simplifications. 

In  our  case,  the  simplifications  concern  the  yaw  of  repose.  A  close  study  of  Eq.(4)  reveals 
that  Qe  is  an  impossibly  intricate  function  of  six  of  the  seven  dependent  variables  (all  but 
A’3)  and  all  five  parameters  C,.  The  magnitude  [oel  is  usually  small  and  only  its  square 
appears  explicitly  (in  Eq.(3),  where  it  constitutes  a  relatively  minor  addendum  to  the  drag 
coefficient).  We  kept  [cve]^  in  our  basic  fitting  equations,  of  course,  but  we  decided  -  with 
a  fairly  clear  conscience  -  to  ignore  its  partial  derivatives  with  respect  to  the  dependent 
variables  and  the  parameters.  It  just  wasn’t  worth  the  formidable  effort  involved  in  adding 
those  complexities  to  already  labyrinthine  expressions. 

The  second  factor  that  came  to  our  aid  was  the  availability  of  an  automatic  derivative- 
taker:  MACSYMA.  This  software  package  is  discussed  in  the  next  section. 


VI.  MACSYMA 

MACSYMA  (project  M.^C’s  SYmbolic  MAnipulation)  is  an  interactive  expert  system 
-  written  in  LISP  -  for  performing  symbolic  and  numerical  operations.  As  elementary 
examples  (from  Ref.  3): 

•  if  asked  for  di.'^in  x)/dx,  MACSYMA  answers:  cosfr), 

•  if  asked  for  f.sin(x)  dx.  MACSYMA  answers:  -cos(x), 

•  if  asked  to  'simplify'  sirv’x,  MACSYMA  answers: 

sin{bx)  5sm(3x)  ^  sin(x) 

16  16  ’’’  S 

M.4CSYMA  can  expand  a  given  function  in  a  Taylor  or  Laurent  series,  invert  a  matrix 
of  symbols,  solve  a  system  of  nonlinear  algebraic  equations,  manipulate  tensors,  solve 
differential  equations,  etc.,  etc.:  the  list  goes  on.  A  very  useful  program.  In  particular. 
M.4CSYMA  was  used  by  us 
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•  to  fiiul  expressions  for  the  partial  derivatives  of  Fj  with  respect  to  the  seven  dependent 
variables  and  live  fitting  parameters,  as  required  in  E(i.(32). 

•  to  convert  those  expressions  into  FORTRAN  statements  for  inserticjn  in  the  subroutiiie 
reciuired  by  FINLIE. 

Construction  of  this  subroutine  is  the  only  really  challenging  part  of  using  FINLIE  and 
forming  the  partial  derivative  expre'ssions  is  the  only  challenging  part  of  constructing  the 
subroutine.  Thus  we  see  that  MACSVMA  played  a  key  role  in  the  coding  of  our  program. 

We  will  not  attempt  a  short  course  in  the  use  of  MACS^AIA  (see.  instead.  Refs  2  and 
3).  Two  points,  however,  are  worthy  of  note. 

Firstly,  it  would  be  helpfid  if  MACSVMA  coidd  -  on  its  own  -  rc'cognize  recurrimi 
expressions  and  take  appropriate  short-cuts.  We  have  fovmd  that  its  abilities  in  this  is'gard 
are  limited  (this  may  be  due  more  to  our  inexperience  with  MACS\  MA  than  to  any 
inhert'iit  shortcomings  in  the  program).  However,  with  a  little  nudging,  those  short-cuts 
can  be  imposed  on  MACSVMA.  Consider,  for  example,  that  ubiciuitous  \': 

V  =  -  ir, )2  +  (U2  -  nw  +  (C3  -  U3)- 

=  /( !h  •  U'i  •  !Jr> .  !J6 .  Cl ,  C> .  C.3 ,  C'l ) 

We  see  that  the  partials  of  V  with  respect  to  four  of  the  dependt'iit  varial)les  and  four  of 
the  fitting  parameters  tire  needed: 

d\  _  [( ~  H  1  )C'i  +  ( i/t>  —  II :}  )C  3] 

dij2  V 

and  so  on.  We  have  found  that  if  M.4CSVMA  is  not  told  differently,  it  will  re-deri\'e  such 
expressions  every  time  V'  is  encountered.  The  resulting  FORTRAN  expressions  for  the  par¬ 
tials  of  the  Fj's  could  fill  pages  of  print-out.  To  prevent  this,  we  introduced  dummy  symbols 
for  these  repeating  expressions  (say,  DVV2  for  dV/dij')  and  instructed  MACSVM.A.  to  us(' 
the  chain  rule  where  applicable: 

and  so  cjn.  This  reduced  the  size  of  the  subroutine  considerably  :md  because'  the  suii- 
routine  is  called  often  by  FINLIE  in  its  integration  routine  -  it  reducf'fi  the'  run  time  as 
well. 

The  second  point  of  inte.’rest  is  the  possibility  of  bugs  in  the'  M.A.CSVMA  program. 
.According  to  Rand  (Ref.  3),  ".Any  system  as  large  and  complicate'd  as  M.ACSA'M.A  is 
bound  to  have  sf)me  bugs  in  it.”  After  considering  the*  alternatives,  we  decide'd  to  trust 
.MACSVMA. 


VII,  Test  Projectile 

To  ge'iK'i'tite  traje'ctorii's.  we  hael  to  postidate  some  proji'ctik'.  re'al  or  hypothetical, 
with  its  associate'fl  physical  and  ae'rodynamic  propc'rtie's.  We  chos('  a  hypothetical  projci'tile 


Table  2.  Macli  Dependency  of  the  Test  Case  Aerodynamic  Coefficients 


Mach 

C'do 

Mach 

Cl. 

Mach 

Ca,. 

Mach 

0.000 

.1394 

0.00 

1.6568 

0.00 

4.197 

0.0 

-.0178 

O.SOO 

.1394 

1.11 

1.6568 

0.40 

4.271 

0.5 

-.0145 

0.S50 

.1425 

1.60 

1.7514 

0.80 

4.663 

1.0 

-.0123 

0.S75 

.1507 

3.00 

2.0584 

0.90 

5.039 

2.0 

-.0096 

0.900 

.1693 

0.93 

5.364 

3.0 

-.0079 

1.050 

.3685 

1.00 

4.850 

1.075 

.3871 

1.30 

4.570 

l.lOU 

.3954 

1.60 

4.478 

1.150 

.3943 

3.00 

4.478 

1.300 

.3716 

1.500 

.3458 

1.700 

.3241 

2.100 

.2849 

3.000 

.2106 

(bearing  some  resemblance  to  the  M4S3A1  artillery  projectile)  with  the  following  i)h\>ic:d 
properties: 

i  (diameter)  =  155  mm 
rn  (ntass)  =  46.947  kg 

Ix  =  0.1595  kg  —  m~ 

Of  the  six  aerodynamic  coefficients  involved  in  our  equations  of  motion: 

C'do-  Cd2-  C.Vp,,.  and 

four  were  selected  to  be  functions  of  Mach  number,  as  indicated  in  Table  2.  This  tttlde 
shows  pairs  of  values:  Mach  number  and  corresponding  coefficient  value'.  Our  code  per¬ 
forms  straight-line  interpolation  for  Mach  values  between  two  entries.  The  remr.iuing  two 
aerodynamic  coefficients  were  assumed  to  b<*  constant: 

Cd2  =  4.171 
Cv,..  =  -l.S 

Note  that  exactly  the  same  aerodynamic  behavior  was  assumed  in  the  fitting  eepiations 
as  was  used  in  the  equations  for  generating  the  trajectories.  It  is  unlikc'ly  that  an  actual 
flight  would  duplicat<'  so  jjrecisely  a  projectile's  assumed  data  bas('  of  ai'rodynamics. 


VIII.  Test  Conditions 

For  this  preliminary  study,  we  chos('  a  rather  <'lement:iry  situation:  one  in  which 
the  entire  trajectory  could  reasonably  be  contaim'd  in  a  single  time  interval.  To  iiiaki-  a 
single-interval  fit  feasible,  we  impo.sed  two  conditions  on  the  (’([nations: 


(a)  the  wind  coiupoiunits  H’l  and  II*}  were  conveniently  cliosen  to  he  liinaar  wirh  A’.* 
ov('r  the  entire  trajectory: 

\Vi  =  -6  +  0.002X2  rn/s 

IV2  ^  0 

U3  =  3  +  0.004A'2  m/s, 
where  one  in/s  is  a  little  less  than  two  knots; 

(b)  the  initial  velocity  (659  in/s)  and  (piadrant  elevation  (SOO  mils  =  45  deg)  were 
chosen  so  that  the  traj('ctory  reinaiiu'd  in  altitude  zone  1. 

Tims  we  see  from  Tahh'  1  that  the  temi>eratnre  was  generated  as 

T  —  2SS.15  —  C\H  dt(j  K 

where  C-,  =  —6.5  deg  K/km,  and  the  air  density  was  g«merated  as 


p=  1.225  h 


V2SS.15 


ky  /ni^ 


Of  course',  the  same  initial  atmositheric  value's  (2SS.15  de'g  K  anel  1.225  kg/iir’’)  were-  input 
te)  the  fitting  e'epiations. 

The  latitude  L  for  otir  test  run  was  arbitrarily  s«.'t  at  45  deg  Xorth:  the  azimuth  AZ 
was  zero:  go  was  9.S0665  m/s^.  The  initial  spin  was  taken  te>  l)e  1.335. 3S3  rael/s. 


IX.  Test  Results 

Figures  1.  2.  and  3  she)W  the*  g('ne*r;)t('d  A*i.  A’2.  anel  A’j  value's.  re*spe'ctive'ly.  ve'rsus 
time.  FINLIE's  task  was  to  recove'r  fre)ni  the)se  pre)saicdoe)kmg  e.mrves  the*  atnie)spherie' 
conditions:  air  density,  temperature'  and  wind  history  ewer  the  course  e)f  the'  traje'cre)ry. 
(In  Figure  3.  the  late'ral  me)tie)n  for  zero  wind  is  ;dso  shown:  in  Figure's  1  anel  2.  the  ne)  winel 
e'urves  we)ulel  lie'  almeist  e)n  top  e)f  the'  given  curve's.) 

Tal)le  3  lists  the-  values  ejf  thre-e*  initial  e'enielitienis  anel  five'  parame'te'is: 

•  the'  column  heaeleel  "Ge'ii.  valuer"  lists  the  vahms  usetd  te)  generate  the)se'  A'l.  A'j.  anel 
A:j  value's  pleette-el  in  Figure's  1  te)  3: 

•  the'  e-e)hmm  he'aele'el  "Initial  e'st."  lists  the'  input  e'stiiiic'ite's  use'el  by  FIXLIE  to  start 
the'  fitting  pre)e-e's.s: 

•  the'  e'e)lumn  he'aele-el  "Final  fitte-el  value'"  lists  the*  value's  e)btain('el  by  FI.XLIE  in  the' 
final  ite'ratie)!!  e)f  its  fitting  pre)e-ess. 


Table  3.  Tost  Case  Results 


IC  and  Parameters 

Gen.  value 

Initial  est. 

Final  fitted  value 

A'lo  (m/s) 

405.983 

403.090 

405.975 

A'io  (m/s) 

405.983 

402.745 

405.983 

AXo  (m/s) 

0 

0.030 

0.003 

C\  (m/s) 

-0 

4 

-5.950 

Cb  (1/s) 

0.002 

0.001 

0.00199 

Ci  (m/s) 

3 

-4 

2.982 

C,  (1/s) 

0.004 

0.001 

0.00400 

C5  (deg  K/km) 

-0.5 

-0.05 

-0.507 

(FIXLIE  can  be  instructed  to  fix  the  valu<*  of  any  initial  condition  or  parainet('r  at  its 
input  value.  This  was  done  for  four  of  the  seven  initial  conditions:  the  position  vector  = 
(0,0.0)  and  the  initial  spin.  Ou  =  1335. GS3  rad/s.  We  decided  to  let  FIXLIE  deteriniiu'  the 
initial  velocity  components;  the  first  estimates  shown  in  Table  3  were  obtained  by  haviii”,' 
the  program  take  position  and  time  differences  for  the  first  two  data  points.) 

We  see  from  Table  3  that  the  final,  fittt'd  values  are  in  vfu'v  close  agrei'iueiit  with 
the  values  used  to  gem.'rate  the  traj«tctory.  Xote  that  the  initial  piiramett.-r  estimates 
are  deliberately  poor;  convergence  of  the  fitting  process  didn’t  dt.'pend  on  starting  chjse 
to  the  answer.  Plots  of  the  generated  and  fitted  air  densities  [pj^n  ttnd  pja)  would  be 
indistinguishable  to  the  eye.  Hence  we  choose  a  more  informative  means  of  comptu'istm: 
the  number  of  significant  digits  in  pju  compared  to  pjen-  Let  E^,  be  the  relative  density 
error: 

E/)  =  Kpjen  ~  Pfit)/ Pgen\ 

A  'well-known'  theorc'in  in  numerical  analysis  states  that 


A  number  is  correct  to  at  Ita.ft  .1  significant  digits  when  its  relative  error  is  nor 
greater  than  0.5  x  10  ' . 


To  determine  J  for  our  density  fit,  we  set 

E,  =  0.5  *  10'-^ 


and  solve: 


.7  =  lofjio 


i  30  i 


.4  plot  of  .7  v.i.  time  is  given  in  Fig.  4.  (Of  course,  we  could  have  plotted  the  integer  part 
of  the  .7  <'xpr(‘ssion.  but  a  smooth  (mrve  seemerl  more  attractive.)  The  plot  indicates  that 
the  fitted  density  is  at  its  most  pn'cist'  near  the  start  and  end  of  the  flight  (that  is.  at  the 
low('r  altitudes,  as  w<’  might  expect).  Imt  even  at  the  higlu'r  rdtitudes  there  is  never  k’ss 
than  about  four-digit  agreement. 


The  tempertiture  match  is  more  etisily  ('xpressed.  W('  see  from  T.ible  3  rh.ar  for  this 
test  run. 

T  - 


Tfa  =  -0.5/7  -  (-0.507/7)  =  0.007// 


ilf  <jK 


where  H  is  in  kilometers.  Thus  at  an  altitude  H  of.  say.  4  km.  the  tempc-rature  error  is 
only  0.028  out  of  (288.15  -  6.5  *  4  =  262.15)  deg  K. 

X.  Conclusions 

The  close  ht  that  we  obtained  seems  to  us  (perhaps  in  our  naivete)  a  rt'inarkable 
result.  The  generated  position  values  (Figs.  1-3)  obviously  contained  within  themselves  a 
core  of  data  sufficient  to  allow  an  excellent  reconstruction  of  tlu'  air  density,  sound  speed 
and  winds  along  the  trajectory. 

Of  course,  the  authors  mustn't  liecome  .so  captivated  loy  their  results  that  they  over¬ 
look  the  greatly  simjilified  conditions  iiiuler  which  those  results  were  obtained.  Much  work 
remains  to  lie  done.  For  <'xample,  futun'  extensions  of  this  report  might  consider: 

•  trajectori(’s  that  cover  mor('  than  one  altitude  zone,  so  that  the  generated  density  aiul 
temperature  can  take  take  different  forms  at  different  times  in  the  same  trajectory: 

•  a  series  of  time  intervals,  where  tlu'  density  and  temp<“rature  are  known  only  for  the 
first  interval: 

•  winds  that  are  iionliiu'ar  with  altitude  within  each  time  interval: 

•  random  noise  added  to  the  data:  to  the  position  values,  to  the  af.'rodynamic  coefficient 
values,  etc.; 

•  different  givcui  trajectory  data  (say.  slant  range  and  slant  velocity,  rather  than  A  i .  Ab.  A' 

•  the  complete  six  degrees-of- freedom  <.'(juations  of  motion,  rather  tha.n  the  modifii'd 
point-mass  etjuations. 

The  possibilities  for  making  conditions  more  difficult  for  our  program  ar(.‘  ;dl  too 
numerous.  Still,  we  have  made  a  start. 
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List  of  Symbols 

(gC'D)(V).  Il/Sl 

azimuth  of  the  1-axis,  measured  clockwise  from  North 

k-J  (IS)  (^)’ , 

Coriolis  acceleration,  — 2ai  X  U 

drag  coefficient:  \drag  /orce|  =  (pV^S/2)Cd 

zero-yaw  and  yaw-drag  coefficients: 

Cd  —  Cdo  +  ^021^61^ 


roll  damping  moment  coefficient: 

|ro//  damping  moment]  —  ±(/9V'^S£/2)f(p^/V')Op 


lift  force  coefficient: 

\lift  /orce|  =  /  _  M/.reiCL,, 


static  moment  ^'oefficient: 

\stotic  moment]  =  E[pV'^S(./2)]ae]CM^ 


Magnus  force  coefficient: 

]Magnus  force]  —  ±(pV'^S/2)(c>^/V’)|cre|Cvy,^ 

fitting  parameters  defining  the  wind,  Ecis.(23)  and  (25),  and 
the  temperature  gradient,  Eq.(26) 


/la(G«r) 


iVs] 


dijj/Oqk-  influence  coefficients 


gun  elevation 


relative  error  in  ‘he  fitted  air  density 
RHS  of  the  fitting  equations  (30),  j  =  1,2,. ..7 


gravity  acceh'ration 


1^1  at  sea  level;  the  STAT  value  at  45  deg.  N.  latitude  is  9.S()GC5  m/s 
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G 
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Gi.  G'2,  G:i 


h 


(1 


hi 


hM 

H 

H, 
Hb 

I. 
J 

kl 

e 

L 


ni 


Q 

R 

S 


fiat-earth  svst('ni  ooniponcnts  of 

~  -OoGXJR,  1-2X2//?,  0) 


g  +  C ,  gravity  plus  Coriolis  acceleration 


1  +  ^a 


n  1 


m/s^] 


flat-earth  system  components  of  G 


l-h^ 


M 


>^1  (liff)  (jf) 


geopotential  altitude.  Eq.(lS) 

value  of  H  at  the  starting  time  of  a  fittin  interval 

value  of  H  at  the  bottom  of  a  STAT  altitude  zone  (Table  1) 

axial  moment  of  inertia 

number  of  significant  digits  in  the  fitted  air  densdy 

R/mi^ 

reference  length 

latitude  at  the  launch  point  (for  Southern  Hemisphere 
firings,  replace  I  by  -L) 

projectile  mass 

the  seven  initial  conditions  (k  =  1,2,. ..7)  and  five 
parameters  (k  =  S,9,...12)  of  the  fitting  equations  (30)  • 

constant  in  the  STAT  air  density  formulas,  (21)  and  (22) 

effective  radius  of  the  earth  (6  356  766  m) 

reference  area,  /4 
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19G2  U.S.  Standcircl  Atmosphere,  Ref.  6 
time 

temperature,  degrees  kelvin 

temperature  gradient  dT/dH;  fitting  parameter  C5 

t(?mperature  at  the  starting  time  of  a  fitting  interval.  Eqs.(27)  and  (29) 
STAT  temperature  at  altitude  Hg,  see  Table  1 
projectile  velocity  with  respect  to  the  earth 
fiat-earth  system  components  of  U 

\v\ 

U  —  IT.  projectile  velocity  with  respect  to  the  air 
speed  of  sound 

flat-earth  system  components  of  V 

w'ind  velocity  with  respect  to  the  earth 

flat-earth  system  components  of  W 

projectile  position  with  respect  to  the  earth 

flat-earth  system  components  of  A’ 

A'l  ;  down- range 

X2  ■  the  height  above  sea  level 

A:j  :  lateral,  positive  to  the  righi  looking  down  rang*' 
dependent  variables  of  the  fitting  equations  (30),  j  =  1.2... .7 
geometric  altitude.  Eq.(16) 


Or 

P 


the  yaw  of  repo.se,  Eq.(4) 
air  density 


p.l 


PB 

0 


'^'2t  ^’3 


air  density  at  the  starting  time  of  a  fitting  interval.  Eq.(29) 
STAT  air  density  at  Hb,  see  Table  1 
axial  spin  rate 
angular  velocity  of  the  earth 
1^1 

fiat-earth  components  of  uJ: 

■jJe{cosL  cos  AZ,  sinL,  —cosL  sin.dZ) 


d{  )/dt 
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